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ABSTRACT 



During the late stage of planet formation when Mars-sized cores appear, inter- 
actions among planetary cores can excite their orbital eccentricities, accelerate 
their mergings and thus sculpture their final orbital architecture. This study 
contributes to the final assembhng of planetary systems with N-body simula- 
tions, including the type I or II migrations of planets, gas accretion of massive 
cores in a viscous disk. Statistics on the final distributions of planetary masses, 
semimajor axes and eccentricities are derived, which are comparable to those 
of the observed systems. Our simulations predict some new orbital signatures 
of planetary systems around solar mass stars: 36% of the survival planets are 
giant planets (> lOM^). Most of the massive giant planets (> SOMq) locate at 
1-lOAU. Terrestrial planets distribute more or less evenly at < 1 — 2 AU. Planets 
in inner orbits may accumulate at the inner edges of either the protostellar disk 
(3-5 days) or its MRI dead zone (30-50 days). There is a planet desert in the 
mass-eccecntricity diagram, i.e., lack of planets with masses 0.005 — O.OSMj in 
highly eccentric orbits (e > 0.3 — 0.4). The average eccentricity (~ 0.15) of the 
giant planets (> lOM^) is bigger than that (~ 0.05) of the terrestrial planets 
(< IOMq). a planetary system with more planets tends to have smaller planet 
masses and orbital eccentricities on average. 

Subject headings: Methods: N-body simulations- planetary systems: formation- 
planetary systems: protoplanetary disk 
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1. Introduction 

To date, around 490 exoplanets have been detectecjl], mostly by Doppler radial velocity 
measurements(e.g., see Udry & Santos 2007 for a review of their statistics). The masses of 
the planets range from order of Earth mass (M^) to tens of Jupiter masses (Mj). Among 
them, the most recently observed HD 10180 system records the most numerous planets in 
exoplanetary systems(Lovis et al. 2010). A study of multi-planet systems is helpful for 
understanding the formation history of planetary architecture due to planetary interactions 
(e.g., Wittenmyer et al. 2009). Recent statistics of single and multiple planetary systems 
have revealed several new possible signatures(Wright et al. 2009): 

1. Including systems with long-term radial velocity trends, at least 28% of known 
systems appear to contain multiple planets. Thus multi-planet systems seem to be 
common. 

2. The distribution of orbital distances of planets in multi-planet systems and single 
planets are inconsistent: single-planet systems show a pileup at period of ~ 3 days 
and a jump near 1 AU, while multi-planet systems show a more uniform distribution 
in log-period. 

3. Planets in multi-planet systems have somewhat smaller eccentricities than single 
planets. 

4. Exoplanets with their minimum masses bigger than Jupiter have eccentricities broadly 
distributed across < e < 0.5, while lower mass exoplanets exhibit a distribution 
peaked near e = 0. 



http:/ /exoplanet.eu 
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5. There may be a positive correlation between stellar masses and the occurrence rate of 
Jovian planets within 2.5AU (Johnson et al. 2007). This might be a reflection of the 
correlation of stellar masses with their circumstellar disks. 

Concurrently, our theoretical insights into the process of planet formation have 
improved greatly. According to the conventional core accretion scenario, planets formed in 
circumstellar disks. Through sedimentation of dust and cohesive collisions of planetesimals. 
Mars-sized embryos will form through runaway and oligarchic growth (Safronov 1969, 
Kokubo & Ida 1998). Angular momentum exchanges between the embryos and the gas disk 
will cause a fast inward migration (type 1) of the embryos in an isothermal disk(Goldreich 
& Tremaine 1979; Ward 1997; Tanaka et al. 2002). In some circumstance, for example, in 
a region that magnetorotation instability (MRI) is active (Laughlin et al., 2004, Nelson 
& Papaloizou 2004), or in a radiative disk (Paardekooper & Mellema 2006, Kley et al. 
2009), the speed of type I migration can be greatly reduced or even with its direction 
being reversed. Thus the embryos can be effectively sustained. As long as they grow 
beyond some critical masses (~ lOM®), quasi- hydrostatic gas sedimentation begins in the 
Kelvin- Helmholtz timescale (Pollack et al. 1996, Ikoma et al. 2000), which may last for 
million years until the final runaway gas accretion sets in. Type I migration of planetary 
cores is helpful to shorten the timescale of giant planet formation, provided a factor of ~ 10 
reducing of the migration speed (Ahbert et al. 2005). The newly formed giant planets then 
undergo type II migration before the gas disk is depleted. 

According to the above scenario of single planet formation, planetary population 
synthesis has reproduced successfully most of the statistical characteristics for the observed 
systems. In a serious of works, Ida & Lin (2004a, 2004b, 2005 and 2008) predicted a 
planet desert with masses in 10-100 M®, due to the runaway gas accretion of giant planets, 
and confirmed the observed correlation that the planet detection rate increases with the 
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metallically of their host star (Fischer & Valenti 2005). Other studies with population 
synthesis of planets, such as Kornet & Wolf (2006), Mordasini et al. (2009a,2009b), give 
more constrains of planet formation model based on observational data. 

However, planets formed in an environment with many siblings. Perturbations from 
neighbor protoplanets may excite their eccentricities, which is helpful to planet formation 
in some cases. For example, the inward type II migration of gas giants in outside orbits 
may trap the planetary cores in inner mean motion resonances (mainly 2:1 resonance), thus 
results in their mergings into hot super-Earths (Zhou et al. 2005, Raymond et al. 2006). 
While planetary perturbations may inhibit the core growth in other cases, e.g., the lack of 
planets in the location of asteroid belt may be due to Jupiter perturbations. Moreover, 
planet-planet scattering is thought as the major cause for the eccentricities of the observed 
planetary systems (Rasio & Ford 1996, Zhou et al. 2007, Chatterjee et al. 2008, Juric & 
Tremaine 2008). Planet population synthesis method without including the planet- planet 
scattering procedure is incomplete to account for the observed planetary signatures. 

By taking the planet-planet scattering effect into consideration, N-body integration of 
the full planetary equations is necessary. Based on N-body simulations and a self-consistent 
disk evolution, Thommes et al. (2008) investigated the formation of giant planets in the 
context of disk evolution, and revealed some interesting tendencies that relate the planetary 
systems with their birth disks. Ogihara & Ida (2009) studied the formation and distribution 
of terrestrial planets around M dwarf stars, and find the final configurations of planets 
depend on the speed of type I migration. 

In this work, we employ the N-body method to study the final assembly of planetary 
systems. The aims of this study are twofold: (1) to explain the observed statistics of planet 
masses and orbit parameters, thus to constrain the parameters that are most suitable for 
planet formation; (2) to guide the future observation by extending our study to the mass 
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regime that is not yet detectable. We adopt the standard model of core accretion scenario, 
which includes the standard type I and II migrations of planets. 

For N-body simulations of planetary formation and evolution, one of the major 
difficulties is the initial conditions of planetary embryos and the embedded protoplanetary 
disk. At the later stage of planet formation, most of the embryos are assumed to be 
formed with their masses ranging from Mars to several Earths (or even bigger). However, 
their radial distributions are quite uncertain. Also, the parameters of viscous disks (their 
extensions, surface density distributions, depletion timescales, etc.) are quite uncertain. 
These parameters may be coupled, and their effects are different. As shown in Thommes 
et al. (2008), disk properties play a key role in determining the planetary migration. They 
chose the mass and the viscosity of the protostellar disk as two free parameters. To simphfy 
the model used in this paper, through some tests, we fix most of the parameters to some 
fiducial values, except letting the gas disk depletion timescale {Tdisk) a-^d disk mass (Eg) 
as two free parameters. The effects of other parameters on the final architecture of the 
planetary systems are also discussed. 

The arrangement of this paper is as follows. First we describe the model and method 
of the paper in §2,. Then some analytical estimation about the orbital configuration under 
type I and II migrations is given in §3. The effects of disk parameters on the evolution and 
final architecture of the planetary systems are discussed in §4. In §5 , the distributions of 
planet parameters (masses, semimajor axes, eccentricities) from simulations are compared 
with observations. To see if our results are credible in two-dimensional model , we survey 
the differences in three-dimensional model in §6. The conclusions are presented in §7, with 
discussions and some implications for future works. 
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2. Model and Initial Setup 

Our investigation starts at the stage that most of the planetary embryos have cleared 
up their feeding zones and obtained their isolation masses. This may correspond to an 
epoch of ~ 1 Myr after the birth of the protostar. Embryos in inner obits might undergo 
type I migration, which is assumed to be stalled at some locations such as the inner edge of 
MRI dead zone of the gas disk (Kretke & Lin 2007, Kretke et al. 2009). For embryos outside 
the snow line, some of them are large enough (several Mq) for the onset of efficient gas 
accretion. They begin to accrete gas, open gaps, and undergo type II migration. Detailed 
model and parameters that we used in this paper are presented as follows. 



2.1. Viscous Disks 

Pre-main-sequence stars accrete gas through circumstellar disks. Although the origin 
of the disk viscosity is still in controversial, the onset of MRI helps to the transportation of 
angular momentum through the disk (Balbus & Hawley 1991). In this paper, we adopt the 
ad hoc a-prescription (Shakura & Sunyaev 1973) to model the effect of mass transportation 
during the stage of classical T-Tauri stars. In such a disk, the kinematical viscosity is 
expressed as = acgh, where Cg and h are the sound speed in mid-plane disk and the 
density scale height, respectively (Table 1). The surface density at stellar distance a is given 
as (e.g., Pringle 1981), 



where M is the gas accretion rate and is the stellar radius. For a T-Tauri disk with 
stellar mass of O.2M0 < < 2.OM0, the observed accretion rate is (Natta et al. 2006, 
Vorobyov & Basu 2009) 

M = 2.5 X 10-«Moyr-^ ij^j . (2) 
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In an optically thin disk, adopting the parameters in Table 1 and let a ^ i?*, we get Eq.(|3]) 
from equation ([T]), 



S, = 2400 gcm-2 -— — - . (3) 



'0. 

For an a-disk with constant a = 10^^, we derive the minimum mass solar nebular 
(MMSN) except with a different slope, (3 = - In Eg/ In a (c.f., /3 = 3/2 in MMSN, Hayashi 
1981, Ida & Lin 2004a). We also assume that the protoplanetary disk has a layered structure 
as a sandwich: inside a location a^rit, the protostellar disk is thermally ionized partly, while 
outside Ocrit only the surface layer is ionized by stellar X-rays and diffuse cosmic rays, 
leaving the central part of the disk a highly neutral and inactive "deadzone" (Gammie 1996). 
The viscosity in the MRI active region could be one or two magnitudes larger than that of 
the dead zone(Sano et al. 2000). By assuming a constant accretion rate across the disk at a 
specific epoch in equation (1), a positive density gradient is expected near Ccrit, which helps 
to halt the embryos under type I migration(Kretke & Lin, 2007, Kretke et al. 2009). 

To model this effect, we let omri and adcad denote the a- values of the MRI active and 
dead regions, respectively. The effective a for the disk is modeled as (Kretke & Lin, 2007): 

"eff(a) = ^ i^^^TTl ) + M+ "MRI, (4) 

where erf is the error function, O.lOcrit is thought as the width of the transition region. In 
this paper, we adopt omri = 0.02,adead = 10~^ (Sano et al. 2000). The small value of 
ttdead is adopted to obtain a reasonable timescale of type II migration (see Eq.[T6] later). 
Alternative choices of Odead are also discussed in §4.1. 

The location of the inner edge of the MRI dead zone, Ocrit, varies with the disk 
temperature, kinematics and mass accretion rate, etc. Here we adopt the expression from 
Kretke et al. (2009), 

„ _ nifi ATT/ ^ N4/9/^*Nl/3/ «MRK -i/5 ..x 
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During the evolution of a T-Tauri star and its disk, M decreases from < 10~^ to ~ 10~^ 
Mq yr~^, according to the infrared excess observation(e.g. Gullbring et al. 1998). So the 
location of Ocrit is migrating inward as time proceeds. To simplify the procedure, we let the 
gas disk depletes uniformly in a timescale Tdisk, so M = Mq exp(— t/rdisk)- Assuming a disk 
mass of O.O2M0 (Andrews & Wilhams 2005), we obtain Mq = O.O2M0/rdisk. 

Substituting the previous equations ^ and ([5]) into equation (j3]), we obtain the 
following surface density for the circumstellar disk: 

where fg is the gas enhancement factor, Sq = 280g cm^^ is adopted in this work so that 
the total mass of the disk up to 100 AU is 0.02Mq for fg = 1, corresponding to the 
average disk mass in Taurus-Auriga star formation region(Beckwith & Sargent 1996). For 
such a disk, the Toomre's criterion for the onset of gravitational instability is expressed 
as Q = = 340(a/lA[/)"3/^/^"^ So Q > 1 holds and gravitational instability will not 
occur up to 100 AU as long as < 10. 

We truncate the inner disk at the disk cavity due to the stellar magnetic field. Around 
the corotation radius, the stellar magnetic torque acts to extract angular momentum from 
the disk and spins down disk material. At the location where the stellar magnetic field 
completely dominates the disk internal stresses, sub-Keplerian rotation leads to a free-fall 
of disk material onto the surface of the star in a funn el flow along magnetic-field lines. 



which results in an inner disk truncation (jKonigl 199ll ). The maximum distance of the disk 



truncation is estimated at ~ 9 stellar radii. Considering that the radius of protostar is 
generally 2-3 times larger than their counterparts in the main sequence stage, the inner disk 
truncation would occur at < O.lAU. Thus we set the inner edge as 0.05AU in this paper. 
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2.2. Isolated Embryos and Type I Migration 

According to the core accretion scenario of planet formation, planetary embryos form 
through cohesive collisions of heavy elements near the midplane of the circumstellar disk. 
Within the solid disk an embryo will grow until it accretes all the dust material round 
its feeding zone so that an isolation body is achieved (Kokubo & Ida 1998,Kokubo & Ida 
2002). At the later stage of planet formation, inward migration of planetesimals under 
gas drag and planetary embryos under type 1 migration may change the distribution of 
embryos, which results in an uncertainty of the initial embryos masses. On the other 
hand, the MMSN model for the sohd disk was derived according to the final solid masses 
of the present planets in the solar system (Hayashi 1981, Ida& Lin 2004a), which reflects 
the final distribution of the embryo, so we adopt the isolation masses with distribution as 
that of MMSN but with an enhancement factor fa — fg- In a gas disk, the instabihty of 
the isolation masses is suppressed until the gas disk is depleted. Assume the dynamical 
instability timescale of a group of isolated masses is the same as the disk dispersal timescale, 
the isolation masses are correlated with their mutual separations, and are given by (Kokubo 
& Ida 2002, Zhou et al. 2007) 

M.,. = 0.16Me ihr,^^^f"(-^f'\^r" (^)''' , (7) 

where is the heavy elements enhancement factor over MMSN model, 77ice = 1 inside 
Oice and 4.2 outside aice, with Oice the location of the snow line, beyond which water is 
condensed as ice from disk gas (T ~ ITO-ft'). fciso is the separation of embryos scaled by 
Rh = (2Miso/3M*)V3a. For Tdisk = 1-10 Myr, Mso = 8 - 10 at lAU and 7 - 9 at 10 AU 
for a disk with fd — i (See Fig. 11 of Zhou et al. 2007 and details therein). 

Due to the variation of the disk accretion rate and stellar radiation at different epoch 
of the protostar, the location of snow hne varies (Garaud & Lin 2007). For a star with mass 
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< 3 Mq, we adopt the snow line location as Kennedy & Kenyon (2008): 



2.7 AU(— )^/9 



© 



IO-8M0 yr- 



t) 



2/9 



(8) 



To model the type I migration of the embryos, we adopt the expression of the migration 
timescale from Cresswell & Nelson (2006), which is also valid for eccentric orbits, 

1 1 fM^\ f fhY 

Tmig,! — 



Ci 2.7 + 1.1^ \MpJ VSgO^ 

1 + (_SZL)5 



(9) 



^ yi.ih) 

where negative (positive) value of Tmig,! corresponds to the inward (outward) migration 
respectively, Mp, r, e, fl are the mass, position, eccentricity, angular velocity of the planet, 
respectively. Due to the MRI effect, /3 = —dlnEg/dlna can be negative near the location 
of maximum pressure Ocrit- C*! is a reduction factor. Lots of literatures have shown that, 
to produce the observed planetary occurrence rate, Ci e [0.03, 0.3] is an appropriate range 
(e.g. Alibert et al. 2005, Ida & Lin 2008). To test its validity in our model, we set 
Ci = 0.03,0.1,0.3 and execute some test runs. As shown in Fig.la, both the choice of 
Ci = 0.03 and 0.1 produce a very slow planet migration embedded in our disk model, which 
may result in the formation of hot Jupiters very difficult. So we set Ci = 0.3 throughout 
this paper. 

The presence of the gas disk will damp the eccentricities of embedded embryos(Goldreich 
& Tremaine 1980). The e-damping timescale can be described as (Cresswell & Nelson 2006), 



''"edap,! 



Qe [M, 



0.78 VM, 



(10) 



where Qe = 0.1 is a normalization factor to fit with the hydrodynamical simulations. 



- 12 - 



2.3. Giant planet formation: gap opening and type II migration 

According to the conventional accretion scenario, giant planets form through three 
major stages(Perri & Cameron 1974; Mizuno 1980; Pollack et al. 1996) : (1) Embryo growth 
stage. Protoplanetary cores form and grow mainly by the bombardment of planetesimals 
before they attain isolation masses. (2) Quasi-hydrostatic sedimentation stage. The 
accretion of planetesimals tapers as their supply in the feeding zone is depleted. This 
induces a quasi-hydrostatic sedimentation and the growth of the gaseous envelope due 
to the loss of entropy. (3) Runaway gas-accretion stage. When the mass of gas envelop 
becomes comparable to that of the core, a runaway stage of gas accretion sets in continually 
until the gas supply is exhausted by either the formation of a tidally induced gap near the 
protoplanet orbit or the depletion of the entire disk. 

Through quasi-static evolutionary simulations, Ikoma et al. (2000) derived the critical 
core mass where significant gas accretion occurs: 

M \ / K ^ 



'10 Weyr J \l cm^ g V 
where Mcore is the rate at which planetesimals are accreted onto the core, and k is the 
grain opacity (see also Rafikov 2006). Due to the uncertainty of Mcore and k, wc assume 
-^core oc fg and adopt Merit = 4M0/°-^^ in this paper. A planet core beyond this critical 
mass will begin to accrete gas in a Kelvin- Helmholtz timescale, tkh — 10^ yr(Mp/M0)~^, so 
that dMp/dt ^ Afp/rxH (Ikoma et al. 2000). However, the accretion rate is also limited by 
the replenishing rate of the materials (Mdisk), so the accretion rate of the planetary embryos 
is expressed as(Ida & Lin 2004), 



— ; — = mm 
dt 



10-9(^)^Meyr-\ Mdisk 



(12) 



When the planet mass grows to sufficiently large, a tidal-induced gap in the gas disk 
forms around its orbit (Lin & Papaloizou 1979). The critical mass for gap-opening is 
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determined by equating the timescale for Type I torques to open a gap (in the absence of 
viscosity) with that for viscous diffusion to fill it in. This gives (Armitage & Rice, 2005) 

With cteff ~ 10~^ in the MRI dead zone, the critical mass is given by: 

Mjjj = 7.5(^)i/'(^)Me. (14) 

After the giant planet opens a gap around the disk, it is embedded in the viscous disk 
and undergos type II migration. As the mass of the planet grows and becomes comparable 
to the disk mass, migration slows down and eventually stops. The migration speed with the 
slow down effect can be expressed as(Alibert et al. 2005): 

da . , 2Ega\ , , 

— = X mm(l, — - — ). (15) 

So the migration timescale is given as: 

= 0.aMy.(|L)-.(^)(^)- 

"^«a;(l,^^^). (16) 

During the migration of the giant planets, their eccentricities will be damped by the 
disk tide(Goldreich & Tremaine, 1979, 1980; Ward 1988), unless the planet is very massive 
(~ 20 Mj, Papaloizou et al. 2001). As the e-damping rate due to the gas disk is quite 
elusive for different mass regimes of the planets, we adopt an empirical formula for the 
e-damping timescale (Lee & Peale 2002) 

Tedap.II = Tinig,Il/K, (17) 

where X is a positive constant with a value ranging 10 — 100. To choose an appropriate 
value of K, we execute some test runs with K — 10, 30, 100. As shown in Fig. lb, the planet 
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eccentricity is damped very quickly with K = 100. For K = 30, the eccentricity can be 
excited and effectively damped at the end of simulation when the gas disk has been nearly 
depleted. Therefor we take = 10 so that some planets in eccentric orbits can finally 
survive according to the observations. 

The gas accretion of a planet will stop when all the gas material around the feeding 
zone(i?H = h) is exhausted, this gives the truncation mass of a gas giant: 

M,,M = m^f\^)M^. (18) 

This limits the mass of a giant planet by pure gas accretion. However, collision and cohesive 
mergings may increase the masses of giant planets up to several Jupiter masses. 



2.4. Equation of Motion 

Based on the model described above, we investigate the late stage formation of the 
planetary system in this paper. The stellar mass is taken as 1 Mq. The gas disk is assumed 
as in Equations ([6]). We further assume that embryos have obtained their isolated masses 
(see Eq.[7]) with fd = fg. We ignore those embryos with mass < 0.1 in the inner and 
outer orbits, so that there are totally 38 ~ 44 isolated embryos initially in circular and 
coplanar orbits in [0.5-13.5]AU for a system with a gas and a solid disk fg = fd = 1- The 
isolation masses and their radial extension will be changed accordingly for different fg = fd- 
The outer boundary of embryos is set according to the rule that, beyond which the core 
masses that can grow within 5 Myrs are less than 0.1 M^, according to the standard core 
growth (Kokubo & Ida 2002, Ida & Lin 2004a). Although this simplification may neglect 
their dynamical frictions to inner massive cores, this effect is similar to the damping effect 
by the gas disk tide, which is more effective and included already in our simulations. 

The angle elements (longitude of periastron, mean motion) of the embryos are 



15 



randomly chosen. Fig. 2 shows an example of the masses and initial locations of the embryos. 
According to equation ([7]), different Tdisk gives slightly different values of fciso and Miso(see 
Zhou et al. 2007). 

The acceleration of an embryo with mass Mj {i = 1, ...N) is given as, 



where and Vj are the position and velocity vectors of Mj relative to the star, the third and 
fourth terms in the r.h.s. of the equation are the accelerations that cause the eccentricity 
damping and migration, respectively. For embryos with Mj < M/ // defined in Eq. ffT^ . 
^mig = ^mig,i (Eq.[n]) and Tedap = Tcdap,i (Eq.[ID]) are used. For those with Mi > Mjji, 
Tmig = Tmig.ii (Eq. [IE]) and Tcdap = Tedap,ii (Eq. [17]) are used instead. Note that, during the 
growth of an embryo, it may cross the critical mass Mjjj, thus pass from type I to type II 
migration, and also the eccentricity-dampping mode is switched. 

In the absence of mutual planet perturbations {Mj = 0, j 7^ z), i.e., without second 
term in the r.h.s. of Eq. ffTOj) . the secular variations of orbital elements for the embryo 
Mj under migration and eccentricity-damping are derived from the classical perturbation 
theory: 




(19) 




TcdapCl + Vl-e^)' 



a 



redap(l + V^T^)' 



(20) 




where uj is the argument of periastron of the embryo orbit. 
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3. Planet Configurations Under Migrations: Analytical Considerations 

Before we present the numerical results, it will be useful to investigate some ideal cases 
that some cores undergo type 1 or 11 migrations only, and to see the configuration of the 
system without considering their mutual perturbations. This will be helpful to understand 
the onset of instability for the full system. 

3.1. Planetary core configurations under type I migration 

Suppose there are N planet cores with masses below M/ // so that they undergo 
type 1 migration. Eq.(|9]) gives aja = r^igj = kM~^aexp{t/Tdisk) for a ^ Ocrit, where 
Mp = Mp/M^, a = a/1 AU and k ^ 0.23Myr if Ci = 0.3 with the expression of from 
Equation ([6]). An embryo with initial semimajor axis Sq will evolve with a = So — ^Mp, 
where ^ = rdisk[l — exp(— t/rdisk)]/^- So the evolution of relative separation between two 
neighboring embryos with mass difference AMp is given as, 

Aa ^ /Aao\ 1 - ^AMp/Aap 
a V «o / 1 - ^Mp/ao ' 

where Afio is the initial separation. For a general case, AM/Aciq = 'yM/ao, and since 
^Mp/cLQ < 1 before a goes to acrit, Eq. ipT]) approximates to, 

--f^)[l + (l-7)eM/«o]. (22) 

So, as long as 7 < 1 , the inward type 1 migration will increase the mutual distance scaled 
by their mutual Hill radii {Rh oc a ). According to Zhou et al. 2007, such an increase will 
enhance the orbital crossing timescale of the system, for example, in the case of isolation 
masses in Eq.(I7j) with 7 = 3/4. However, due to the perturbation of giant planets and their 
mutual perturbations, embryos in inner orbits will have their eccentricities excited, which 
may result in instability. 
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3.2. Giant planet configurations under type II migration 



Similarly to the previous analysis, giant planets undergoing type II migration with 
a timescale of Eq.f lT6|) may modify their mutual separations. With the same notion of 
previous subsection, we assume that the planet mass is much smaller than the inner disk 



k' ^ 0.6 Myr from equation (fT6|) . So a = ao — t/k', and two neighboring planets will have 
constant separation ASq. When it is scaled by Rh oc a at time t, 



This means that their mutual separation scaled by Hill radii will increase during the inward 
type II migration. However, secular perturbations among them will excite their eccentricities 
during their convergent migration, thus destabilize the system if their eccentricities are high 
enough. 



We numerically integrate the equations of planet motion ( JTOl) with a time-symmetric 
Hermite scheme (Kokubo et al. 1998, Aarseth 2003). Regularization technique is used to 
handle the collisions between embryos, all the embryos have their physical radii (a mean 
density of 3 g cm~^ is assumed), and mergings are expected when the mutual distance of 
two embryos is less than the sum of their physical radii. We assume a perfect inelastic 
collision between embryos. The inner boundary of the gas disk is set as 0.05AU. If a planet 
migrates to < 0.04AU, we remove it from further integration. The external boundary of 
the gas disk is set as lOOAU. If a planet evolves to the orbit with a > 50AU, we consider it 
as being escaped from this system. The tidal effects between the host star and the close-in 
planets are not included in our model for their long effective timescales. 



mass, and the disk has a constant acs- Thus one has a/a = r^ig^u = k'a for a ^ ctcrit, where 




(23) 



4. Numerical Results: Dependence on Disk Parameters 
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Figure 2 shows an example of such an evolution. Initially 44 embryos are put, with 
their masses and initial locations shown in Fig. 2 (a). Finally, three planets are left, with 
masses S.IM^, 238. 2M0 and IST.GM^ (from inner to outer). Moreover, the outer two 
planets passed the periastron alignment (with W2 — ^3 librating around 0) during the 
evolution. 

As we mentioned, parameters of the protoplanetary disk are quite uncertain with the 
present ability of observations. So in this section, we mainly investigate the effect of the 
disk mass, viscosity and depletion timescale on the planet evolution. 

4.1. Influences of the disk viscosity (adead) and disk mass (fg) 

The variations of the disk mass and viscosity are modeled by two parameters, i.e., fg in 
Eq. dS]) and adead in Eq.(^. We set fg = 0.3, 1, 3, 10 and adead = 10"^ 10"^ 10"^ to survey 
their contributions to the architecture of planetary systems individually. Five runs for each 
parameter are executed. We fix the surface density profile as well as the disk depletion 
timescale {r^isk = 2Myr) for these test runs. 

In our model ttdead mainly works on type I migration via Eq.dH]) and type II migration 
via Eq.f lT6|l . Since oc (oe//)"^ (Eq. |6]), thus r^igj oc Oe//, i-e., planets in disks 
with smaller aeff{^ 10~^) migrate faster, which results in a smaller average semimajor 
axis(Fig.3a). For type II migration, if the planet mass is small (Mp < 2'Lga?, which is 
about 9OM0 at 5AU and ITOMq at lOAU), the migration speed decreases with Oe//, thus 
larger Uejf (~ 10^^) may also lead to a small average semimajor axis. Thus the disk 
viscosity may change the final location of the planets. So we take a most plausible value 
(ttdead = 10~^) in the following simulations(Sano et al. 2000). On the other hand, the 
average eccentricities are around 0.1 with large variations for all cases, which seems to be 
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insensitive with adead- The reason is that small adead leads to a high density of gas disk, 
resulting in two competing effects: (a) a fast type I migrations so that a large frequency 
of embryo-scatterings is expected, which enhance the eccentricities of embryos, (b) a fast 
e-damping rate with a small Cmean- As the timescales of these two effects are comparable, 
emean IS nearly independent of ttdead- 

There are three ways for fg to influence the orbital evolutions: the surface density of 
the gas disk(Eq. [6]) , initial isolation masses (Eq.[7] with fd = fg), and the critical mass for 
the onset of runaway gas accretion (Eq.[TT]). Larger fg may result in larger cores extended 
to an outer region. As we assume Mcore oc fd, we have M^it oc fg''^^- Comparing with 
the masses of the initial embryos(oc fg^'^), it's much easier to form massive gas giants in a 
massive gas disk. 

As shown in Fig. 3b, when fg > 1, a^ea.n are all around 3 — 4AU. Although the 
migrations should be faster for planets in disks with larger fg, the initial embryos extended 
to an outer region in these cases, which results in similar averaged locations. For a less 
massive disk(/g = 0.3), these small embryos initially in the inner region have never accreted 
gas, and thus they experienced type I migration throughout the evolution, therefore they 
have smaller amean and Cmean- More effective eccentricity-dampping results in a smaller 
emean for /g = 3 than that of fg = 1. When fg = 10, the largest embryos can reach SOM® 
via collisions in IMyr. Finally the planets left in these systems are mainly gas giants, and 
their eccentricities are damped much less effectively than those of small planets. This is 
why Cmean of fg = 10 is much larger than that of fg = 3. 
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4.2. Correlations with the disk Ufetime ( Tdisk) 

The role of the gas disk on the planet evolution is manifold: it causes the inward 
migration (type I or II) of protoplanets as well as the tidal damping of their orbital 
eccentricities. A long survival gas disk may be helpful to the formation of giant planets, 
while planets in a short-lived disk may not have enough time to accrete gas, thus remain 
either terrestrial or Neptunian planets. To investigate the effect of the disk depletion 
timescales (rdisk), we set 11 values of Tdisk, evenly ranged from 0.5 Myr to 5Myr in a 
logarithm scale (Haisch et al. 2001). For each Tdisk we did 20 runs of simulations by choosing 
the orbital phase angles of the embryos randomly. So totally we did 220 simulations for 
fg = I. All the simulations are stopped at i = 10 Myr. 

An interesting problem for planet formation is the growth epoch of different types of 
planets. We define (somewhat arbitrary) the following two types of planets: gas giant planets 
(GPs), includes massive GPs (Mp > 30M^) and Neptune-sized GPs (lOM® < Mp < SOM®); 
terrestrial planets (TPs), including Super-Earth TPs (IM® < Mp < lOM®) and Sub-Earth 
TPs(Mp < IM0). Fig.4 shows the distribution of planet semimajor axes for the 220 runs 
of simulations at some epoches with fg = 1. Only TPs arc present at t = lO'^yr (Fig. 4a). 
When t = IMyr, runway gas accretion occurred, thus a few GPs appear (Fig.4b). The most 
efficient growth of GPs occurred at 1 — 3Myr, thus the number of GPs increased very fast 
(Fig. 4c), until they become the dominant members at the end of simulations(Fig.4d). 

Fig. 5 plots the correlations between the properties of the final systems and Tdisk- 
Basically there are two regimes. For short-lived disks with Tdisk < IMyr, the forming 
planets have small average masses (70 — 200 M®, Fig. 5c), the average semimjor axes of these 
systems are large (5 — 6AU, Fig.5d), indicating most of the surviving planets were formed in 
distant orbits, while migrations have not affected their orbital architectures strongly. The 
average eccentricities (~ 0.15) are relatively big (Fig. 5a), showing the tidal damping effect 
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is not very effective due to tfie sfiort lifetime of the disk. On the contrary, planetary systems 
with a longer disk lifetime {Tdisk > IMyr) tend to have larger average planet masses (~ 1 
Jupiter mass ), with lower averaged eccentricities (~ 0.1) due to the long period of the disk 
damping. The semimajor axes of these systems decrease from 4 AU to 2 AU, indicating 
inward migration indeed plays an important role in sculpting their orbital configurations. 

Systems with Tdisk < l.SMyr do not have much differences on the surviving number of 
the planetary systems, while A^mean decrease slightly as Tdisk > l.SMyr increasing (Fig. 5b). 
The maximum of A^mcan at ~ l.SMyr, is due to the mechanism of halting small planets 
near Ccrit (see §5.1). In disks with smaller Tdisk, the effects due to gas disks are not so 
efficient, and the surviving number of planet systems mainly depends on the interactions 
between embryos. In the disks with longer lifetime, giant planets experience sufficient type 
II migration and will be closer to the boundary of MRI region, where small planets were 
halted (see §5.1). Thus small planets are easier to be scattered out of the system or hit the 
host star, leaving fewer surviving planets. 

One of the major effects that N-body simulations can describe is the eccentricity 
excitation due to mutual planetary perturbations. Such excitation may lead to some 
embryos being scattered out of the system, while others may be merged. Thus the 
final number of the planets should be greatly reduced from that of the initial embryos. 
Correlations between the number of survival planets and the averaged mass as well as 
eccentricity of the planets in a planetary system are shown in Fig. 6. They are fitted by 

emean = 0.65 X 0.67^'««, 
M^ean = 1.28M^ X 0.85^"'ft. (24) 

These correlations show that, in multi-planet systems like solar system, planets basically 
have relative lower averaged masses and eccentricities than those with single or few planets. 
Qualitatively, this law is easy to be understood. In order to achieve a longer orbital crossing 
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time in a multi-planet system, either the eccentricities of the planets must be small, or 
the planets should have large mutual separation scaled by their Hill radii, thus smaller 
planetary masses is helpful to achieve a stable system(Zhou et al. 2007). 

5. Comparing with Observations 

To compare our results with the observations, we set the mass of the disk (fg) and its 
depletion timescale (rdisk) as two free parameters. We take 11 values of Tdisk evenly ranged 
from 0.5 Myr to 5Myr in a logarithm scale(Haisch et al. 2001). For each Tdisk, "we choose 
fg = 0.3, 1, 3, 10, and execute 6,10,5,3 simulations, respectively, to fit for the Gaussian 
distribution of log {Mdisk/ Mq) with a mean value of = —1.66 and a standard deviation 
a = 0.74, according to the observations of Taurus- Auriga (Mordasini et al. 2009a). So 
totally we execute 11 x 24 = 264 runs of simulations. 

We also make some statistical plots from the observed system]^. To show planets that 
may be not observable yet, we distinguish planets of our simulations with detectables (with 
the induced stellar radial velocity V^^ > 3 m s~^) and undetectables (Vr < 3 m s~^). Taking 
the mean value of sini = 0.6, among the 1437 survival planets, 959 planets (66.7% of total 
) from our simulations are undetectable, which contains a large number of small planets 
(< IMe). 

5.1. Semimajor Axis Distributions 

Fig. 7a and b show the semimajor axis distribution from both the observations and the 
simulations. The observed distribution shows two peaks(Fig.7a). (1) At 1-3AU. This is 
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roughly the snow hne (where water is frozen) of the system, where planet cores may be 
stalled under type I migration (Kretke & Lin 2007, Ida & Lin 2008), subsequent accretion 
of gas makes them giant planets. (2) At 0.04-0.06AU (or 3-5 days). This is roughly the 
inner edge of the gas disk, where planets under type II migration will be stopped (Lin et al. 
1996). From our simulations there are lots of planets beyond 3AU, thus the lack of planets 
at > 3 AU is due to observational bias. 

Besides, we show that there is an extra pile-up of planets: (3) at around 0.2AU (or 
around 30 days), which has not been revealed by observation yet. This location corresponds 
to the inner edge of the MRI dead zone (ocrit) where small planets under type I migration 
may be halted. However, as the location of Ocrit (See Eq.[8]) moves inward, in the case 
of short Tdisk, its migration speed is faster than that of the type I migration, therefore it 
is difficult to halt these planets near acrit- When Tdisk ~ l.SMyr, ^ |mig,i~ ^^f^; so if 
Tdisk > 1.5 Myr, the mechanism of halting small planets near a^it works. After r^isk, the 
migration rate is reduced due to gas depletion. In fact, Ocrit ~ 0.2 AU at t = 2rdisk- This 
explains the accumulation of planets near 0.2AU. Some small planets were scattered into 
0.04AU. 

Our simulations also show that, the majority (70%) of giant planets is located in 
1 — lOAU. Most of them are massive GPs and experienced type II migration. Due to the 
depletion of the gas disk, they can not migrate to the proximity of the star, hence they 
halted outside 1 AU. Due to observational bias of radial velocity measurements, only those 
massive giant planets at < 4 — 5AU are revealed. 

5.2. Mass Distributions 

Our simulations show mainly there are two peaks for planet masses: (Fig.7c). 
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(1) At l-2Mj. They are Jupiter-sized giant planets that have grown up with efficient 
gas accretion. Although our simulations reproduced this peak, there are not enough amount 
of massive planets at Mp ~ 5 — 8Mj. 

(2) At 0.1 — 3M0. There are a large number of terrestrial planets, which have grown 
up mainly by mutual collisions under the perturbations of outside giant planets. This peak 
has not been revealed by the observations yet, and can be checked by future observations of 
terrestrial planets. 

One more small peaks are revealed by our simulations: 

(3) Around lOMj. Massive embryos (M ~ SOMq) had been formed after IMyr in 
disks with fg — 10 and only experienced time-limited type II migration. More massive 
giant planets(> lOMj) may be born with some other mechanisms, e.g., the gravitational 
instability scenario. 

The desert at 10 — 2OM0 are Neptune-sized planets. According to our model in Section 
2.3, the Kelvin-Helmholtz timescale are quite short(~ 10"^yr) for this mass regime and 
the runaway gas accretion makes them become giants quickly. Only a few planets in the 
inner region with Miso ~ 10 — 2OM0 survive in this desert. Also merging from lower planet 
embryos is possible to form Neptune-sized planets. 

5.3. Eccentricity Distributions 

The distribution of eccentricities from simulations is similar to that from observations. 
In our simulations, the eccentricities of survival planets vary from to 0.84, while the 
observations show a maximum eccentricity of about 0.9 (Wright et al. 2009). We fit both 
distributions (simulations and observations) with an exponential decay function in the form 
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of P{e)de = N{e)/Ntot ~ exp{—Ae)de. Since P{e)de = 1, we have 

1 — exp[—A) 

where A = 4.2 for observations and A = 7.8 for simulations of planets with Vr < 3m/ s, 
as shown in Fig.7e and Fig.7f. The larger value of A from simulations indicates a steeper 
slope, showing more planets with small eccentricities are still not detected, as shown in 
Fig.Sd later. 



5.4. Correlation Graphes between a, e and Mp 

Figure 8 presents the correlation diagrams of a, e and Mp of our simulations (right), 
with a comparison to the observational data (left). Our simulations reproduced all the 
three correlation plots quite well. Especially there is a gap in the Mp — e plot, indicating a 
planet desert (0.005 ~ 0.08Mj) depends on the eccentricity (Fig. 8c & d). Fig.8d also shows 
a tendency that giant planets {Mp > IOM0) tend to have lager eccentricity on average. To 
show this more clearly, we reproduce the eccentricity-semimajor axis correlation plots and 
the eccentricity distribution plots for giant planets (Mp > lOM^) and terrestrial planets 
{Mp < IOM0) in Fig. 9. Giant planets have average eccentricities ~ 0.15 at all locations 
(except > lOAU), while terrestrial planets have average eccentricities ~ 0.05. Although 
both of the e-distribution can be fitted by exponential law in Eq.f l25p . their coefficients A 
are different : 7.78 for Mp > IOM0 and 21.0 for Mp < lOM©, indicating small mass planets 
tend to have smaller eccentricities. 

To understand the Mp — e desert, i.e., very few planets with masses 0.005 ~ 0.08Mj 
have eccentricities larger than 0.3 — 0.4, we plot the e-damping timescales Tedap as a function 
of planet masses at 0.3, 1, 3 and lOAU in Fig. 10. For a planet, when Mp < Mj ij, Tedap is 
calculated by Eq. (ITU]) using a mean eccentricity e = 0.05 and Sq = 280 g cm~^. Otherwise 
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Tedap is calculated by Eq. f|T7|) with viscosity adead = 10""^. Due to Eq. f lTOj) . small planets 
have longer T^dap- The jumps of r^dap occur at Mp = Mfjf. r^dap keeps horizontal until the 
mass of the planets becomes comparable to the disk mass (Mp > 2X^0^ = 21M0(a/lAU)). 
In the braking phase of type II migration, more massive planet have a longer e-damping 
timescale according to Eq. ( fTTj) . Assuming an effective damping of eccentricity if Tedap 
is less than ~ IMyr, we obtain a mass regime 0.005 ~ O.OSMj which is consist with the 
Mp — e desert from both observations and our simulations in Fig. 8c & d. As shown in 
Fig. 10, Tedap of massive GPs is in the horizontal region or beyond this, so generally they 
have longer e-damping timescales than those of the TPs. Therefore, the mean eccentricity 
of TPs is smaller than that of GPs. 

6. Simulations in Three-Dimensional Model 

In all the above simulations, we adopted coplanar planetary model (2D). When the 
orbital inclinations of the planets (3D) are included in the simulations, their final orbital 
characteristics of the planet systems may be changed due to the different collision timescales 
between 2D and 3D simulations (Chamber, 2001). In this section, we reveal the effects 
of including orbital inclinations with some more simulations. We executed 20 runs in 3D 
model, setting the same initial conditions with those in 2D model except that the initial 
inclinations are set as 1 degree for all embryos, with randomly chosen. The disk parameters 
are also set as a standard value in our simulations: fg = 1, T^isk = 2Myrs and adead = 10"''. 

6.1. Collisions with Different Masses 

During the whole evolution period we simulated (lOMyrs), there are 606 collisions in 
the 20 runs of 2D simulations, while 471 collisions occurred in the corresponding runs of 3D 
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model. In the 2D runs, the number of embryos colhding with the host star or being ejected 
out of the system is 120, comparing with the number of 231 in 3D runs. It is consistent with 
the result of Chambers (2001). As seen in Fig. 11a, the cumulative distribution functions 
(CDF) of collisions in two models show that most collisions in 3D runs occurred in relatively 
later stages. Thus earlier collisions in 2D model produce embryos with larger masses. 
To show the effects of earlier collisions to the final planetary architectures, we divided 
the collisions as three regions, as shown in Fig.llb: Region I (Mq < Merit, Mi < Merit), 
II (Mo < Merit, Ml > Merit) and III (Mo > Merit, Ml > Merit), whcrc Mo represents the 
larger mass of the embryo before a two-body collision, while Mi represents that after the 
collision. Merit = 4M0 is the critical mass beyond which the efficient gas accretion sets in 
(see discussions below Eq.[TTj ). 

Region I: Most collisions (2D: ~ 75%, 3D: ~ 85%) occurred in this region. These 
collisions only influence the Earth-like planets(TPs, see §4.2). Due to the shorter collision 
timescale in 2D model, embryos have larger masses on average during early time and 
experience faster type I migrations and e-damping according to Eq[9] and Eq[TO]. Therefore 
the TPs have a smaller mean eccentricity (c.f., emean = 0.034 for 2D and Cmean = 0.047 
for 3D) and semimajor axis (see Fig. 11c). As a results of more collisions occurred in 2D 
simulations, the final systems are expected to have less number and larger masses of TPs 
than those from 3D simulation. In fact, we find 11 TPs left among the total 74 survived 
planets in the 20 runs of 2D simulations, while it is 24 out of the total 98 survivals in the 
corresponding 3D simulations, with a small average mass (Fig. 11c). 

Region II: About 8.5% (2D) and 8.4% (3D) collisions occurred in this region. These 
collisions makes the embryos with masses < Merit grow beyond the critical mass. However, 
this effect is limited by the small fraction of collisions and the slow gas accretion rate(Eq[T2]) 
for embryos with masses slightly larger than Merit = 4M0. Thus the differences between 2D 
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and 3D models due to collisions in this region can be ignored. 

Region III: Roughly ~ 16.5% (2D) and ~ 6.6% (3D) collisions occurred in this region. 
These embryos can accrete gas before collisions. As type I migrations speed is much faster 
than that of type II, the final locations of the planets are mainly determined by their 
planets cores. As planet cores in 2D model undergo more collisions, they have bigger 
masses thus a fast migration speed, so their final average semimajpr axis is small than 
that in 3D (Fig.lld). As inner region has less gas to accretion, Mgas oc J^gaRn, while 
Rh oc {Mp/3M^y/^a and S oc a ), the final giant planets in 2D simulations have less 
masses(Fig.lld). As the type I eccentricity damping rate is much small than that of 
type II (Fig. 10), and from Eq.( !TO|) . Tedap,! ~ a^/^M"^, thus larger masses of planets in 3D 
simulations have shorter Tedap.i, thus the mean eccentricity (cmean ~ 0.086 by simulations) 
in 3D model is smaller than that (cmean ~ 0.15) in the 2D model. 

6.2. Differences in Statistics 

As most of our results in §5 are in statistics, we focus on the statistical differences 
between 2D and 3D two models. The statistics of semimajor axes and masses of planets in 
two models are presented in Fig. 12. Some peaks and deserts in 2D model are reproduced 
in 3D model, such as the peaks at acru ~ 2AU, M^it and 1 — 2Mj. However, There 
are still some different characteristics, especially the planet deserts at around O.lAU and 
lOMe < Mp < O.lMj in 3D model. The lack of planets around 0.1 AU in 3 D model 
may be due to the prolonged collision timescale, a much longer time of integration may 
results similar qualitative results with 2D model. The deficit of planets with masses 
IOM0 < Mp < O.lMj in 3D model make the planet desert in §5.2 more obvious. 

Fig. 13 shows the orbital inclination distribution for the survival planets in the 20 runs 
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of 3D simulations. As one can see, most of the planets remain in < 2 degrees, only few 
planets with smaller masses have inclinations ~ 10 degrees. 

To summary, we point out that the differences between 2D and 3D simulations are not 
much, most of statistical results in 2D are qualitatively credible. 

7. Conclusions and Discussions 

Based on the standard core accretion model of planet formation, we simulated the final 
assemble of planets by integrating the full N-body equations of motions. The stellar mass 
is fixed at IMq. The circumstellar disk is assumed to have an effective viscosity parameter 
a^s (Eq.[l]), and undergoes an exponential decay in a timescale of rjisk = 0.5 — 5Myr. 
Initially embryos with isolation masses larger than O.IM0 (Eq.[7]) are put in each system. 
For embryos below (or above) the critical mass defined in Eq. f fT^ . they will undergo type I 
(or type 11, resp.) migration. The type I migration of embryos will be stalled near the inner 
edge of the MRI dead zone defined in Eq.([5]). The inner edge of the disk is set as 0.05AU 
where giant planets will be stalled under type II migration. The equations governing the 
embryo's motions are shown in Eq. (JT9|1 . 

We have investigated the influences of different parameters(viscosities, disk masses, 
disk depleted timescales) on the final architectures of the planetary systems. Disk viscosity 
affects the planet locations by determining their migration speeds. In the regime of 
«dead ^ [10^"^, 10^^] , plaucts in disks with small viscosities migrate fast under type I 
planet-disk interactions, which results in a small average semimajor axis. For a moderate 
disk viscosity, the planetary system shows a moderate averaged locations (Fig.3a). The 
average eccentricity does not show obvious correlation with the disk viscosity. Disk mass 
affects the system through the initial core masses. Large planet cores can be formed in the 
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outer region of a massive disk, so that giant planets are easier to form. 

Disk depletion timescale (raisk) plays an important role in final planetary masses and 
orbits. For short-lived disks with raisk < IMyr, the forming planets have small masses 
(70 — 200 M®), with large average semimjor axes (5 — 6AU). This indicates that most of the 
surviving planets are formed in distant orbits, while planet migrations did not affect their 
orbital architectures yet. Due to the short hfetime of the disk, their average eccentricities 
(~ 0.15) are relatively big. Planetary systems with a long-hved disk (rdisk > IMyrs) tend to 
have large planetary masses (~ 1 Jupiter mass ), with low averaged eccentricities (~ 0.1) 
due to the long time disk-damping. The average semimajor axes of these planets are 
ranging from 4 AU to 2 AU, indicating inward migration indeed plays an important role in 
making planets in close-in orbits. 

Comparing the simulations to others, we find the number of planets being trapped in 
mean motion resonances (MMRs) are smaller than those obtained, such as in Terquem & 
Papaloizou (2007). Instead, we get a large amount of planet pairs that have the history of 
passing periastron alignment (See Table 2, with the difference of their longitude periastron 
librating at 0°). The reason that we did not observe many survival planets in MMRs may 
be due to the fact that, most of our survival planets arc giant planets. The choice of a small 
disk viscosity (ctdead = 10~^) makes the type II migration too slow to make a significant 
convergent MMR trapping. For terrestrial planets, the speed of type I migration is faster 
than that of type II migration of the giant planets in outside orbits, so they are not easy to 
be trapped into the MMRs of giant planets as well. 

Statistics of 264 simulated systems reproduced qualitatively the main features on 
the planet masses and orbital parameters for the observed exoplanetary systems. If we 
classify the planets into two major categories: giant planets (CPs), including massive GPs 
(Mp > 30Me) and Neptune-sized GPs (lOM® < Mp < 30Me); terrestrial planets (TPs), 
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including Super-Earth TPs (IM® < Mp < lOMe) and Sub-Earth TPs(Mp < IM®), then 
results from our simulations have the following implications on these planets. 

Occurrence rates of planets. The ratio of GPs relative to TPs is low, i.e., 514 to 923 (or 
36% to 64%) in our total 1437 survival planets (table II). This is mainly due to the fact that 
only those massive embryos can accrete sufficient gas to form GPs. Planets with smaller 
masses are easier to be scattered out by N-body interactions, which is the major difference 
with the population synthesis simulations of single-planet systems. Moreover, there exist 
some correlations between the survival number of planets and the average eccentricity (or 
average planet mass) of a planetary system, i.e., a planetary system with more planets 
tends to have smaller planet masses and orbital eccentricities(Fig.6 and Eq. [21]). These 
correlations are consistent with the stability of the system, i.e., a system with planets in less 
eccentric orbits and with larger mutual separation scalled by their Hill radii tends to have 
longer crossing timescale, thus is more stable (Yoshinaga et al. 1999, Zhou et al. 2007). 

Locations of giant planets (GPs). Most (298, ~ 58% of total 514 ) of GPs locate at 
orbits with semimajor axes 1-10 AU , with only ~ 35% (182) at inner orbits < 1 AU. This 
may be linked with the snow line (~ 2 — 3AU) of the system. Due to the surface density 
enhancement of about 3 — 4 (Eq. [7]), isolated masses beyond the snow line are larger thus 
they are easier to accrete gas and become giant planets. Neptune-sized GPs are formed in 
our simulations also, with its amount being much smaller than massive GPs ( ~ 7.5%, 42 
of 1437), and they mainly locate at two parts of the systems: either at ~ lOAU, which were 
scattered out by first formed giant planets (like the formation of Uranus and Neptune), or 
in the inner edge of the gas disk (0.05AU), which seems to be stalled there under type II 
migration. 

Locations of terrestrial planets (TPs) . TPs are almost evenly distributed, except a 
group lying at 1-2 AU(Fig.4d, for fg = 1), just inside the snow line(~ 3AU). Also there 
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is a slight pileup of planets (both GPs and TPs) at 0.2 — 0.3AU (30-50days), which may 
correspond to the inner boundary of the MRI dead zone. According to Kretke & Lin 
(2007) and Morbidelli et al. (2008), super-Earths are easier to be stalled there from type I 
migration, and grow up to giant planets by gas accretion. 

Eccentricities of planets. There are very few planets with masses 0.005 ~ O.OSMj that 
have eccentricities larger than 0.3 — 0.4, The average eccentricities of giant planets are larger 
(~ 0.15) than those of the terrestrial planets(~ 0.05, Fig.9). According to our simulations, 
the underlying mechanism is the relatively long e-damping timescale of massive planets 
due to the gas disk, especially when the planet mass is larger than the inner disk mass, see 
Eq.(Il6D and Eq.(fT71). 

We compared our results with some 3D simulations in §6. The qualitative differences 
between 2D and 3D models are not big, indicating our conclusions based mainly on 2D 
simulations are still reliable. 

Now let's compare above conclusions to the five new observational signatures that 
stated in §1. According to our simulations, 260 runs in total 264 ones result in multi-planet 
systems, which rates to 98%, comparing with that from the observations (> 28%). Planets 
in multiple systems have smaller eccentricities than single planets, which is revealed clearly 
in Fig. 6. For signature-2, we did not classify our a,e and Mp distributions with single and 
multiple systems, as this classification has some uncertainty concerning undetected planets. 
While Fig.Qa supports signature-3, massive planets seem to have larger eccentricities than 
those of smaller mass planets. The implication of this is interesting. As most of the observed 
exoplanets have Jupiter-sized planets in elliptical orbits, we can expect more planets with 
small masses in near circular orbits of multi-planet system, like terrestrial planets in the 
solar system. Signature-4, i.e., massive planets tend to have larger eccentricities, is revealed 
by our simulations (Fig.9b). For signature-5, as we fix the stellar mass at 1 Mq, it can not 
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be tested in the present work. We will investigate this correlation in future works. 

However, due to some limits and uncertainties of the parameters we chose in this 
model, the conclusions obtained in the paper may be limited. For example, we investigated 
planet formation around 1 Mq stars only, although with different sizes of the sohd and 
gas disks. Not only the disk mass but also the accretion rate as well as metallicity have 
correlation with the mass of the host star. For more broader parameters, the occurrence 
ratio of planets between TPs and GPs may need further investigations. 

One of the major uncertainties arises from the initial masses and distributions of the 
embryos, which are the building blocks of planets. In this paper, we adopt the assumption 
that most of the embryos have already cleared their nearby heavy elements and achieved 
their isolation masses. This assumption is based on the full N-body simulation of planet 
growth (Kokubo & Ida 1998), and will be too ideal when type I migration of embryos are 
taking into considerations. Also, initially there might be some smaller embryos between 
giant planets so that they may be scattered through planet scattering and become planets 
like Uranus and Neptune. 

The second major uncertainty comes from our poor knowledge of type I migration of 
the embryos. The speed (and even the direction) of type I migration will affect the number 
and locations of survival terrestrial planets, especially planets in habitable zones (e.g., 
Ogihara & Ida 2009, Wang & Zhou 2010). We will do more investigations to that end in 
our forthcoming works. 

Thirdly, our simulations mainly focus on the stage that the gas disk is present and 
will be depleted exponentially. Further evolution due to the effect of a planetesimal disk 
was not included. During the later stage of planet formation when the gas disk is totally 
depleted, planet-planetesimal interactions may damp the eccentricities of planets through 
dynamical frictions, and may induce migrations through angular momentum exchanges 
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with embryos and scattered planetesimals (Fernandez & Ip 1984; Malhotra 1993). In the 
solar system, numerical simulations show that, Jupiter will drift inward, while Saturn, 
Uranus and Neptune may migrate outward, resulting in a divergent migration. During the 
migration, the cross of 2:1 MMR between Jupiter and Saturn excites the eccentricities of 
four giant planets (Tsiganis et al. 2005). Such evolution may occur on the timescale of at 
least hundreds of million years, which is beyond the ability of our simulations. 

Furthermore, we ignored the tidal effect between a host star and close-in planets. The 
tidal dissipation of the star-planet system may damp the eccentricity of the planets in 
close-in orbits in a Gyr timescale. Also, the gravitational potential of the gas disk was not 
included in our model. During disk depletion, the sweeping of secular resonance through 
inner region may excite the eccentricities of inner orbits, thus induce further mergings of 
terrestrial planets (Nagasawa et al. 2003). 

Although with many restrictions, the present paper, aiming at the statistics of final 
assembling of planetary systems in the standard formalism (which includes type I and 
type II migration, gas accretion, etc.), reproduces most of the observed orbit signatures. 
Thus we think that the predictions by the simulations are helpful for guiding future planet 
detections. 
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Table 1: Notations used in the paper, M* = M^/Mq,?- — r/lAU, 

T — 280M*f~^/^ disk temperature by stellar radiation 

Cg — 1.2km s~^Ml^^f~^l'^ sound speed in the ideal gas with adiabatic index 7 = 1.4 

Vk — 29.8km s~^Ml^'^f~^^'^ orbital speed of Kepler motion 

h — rcg/vk — 0. 04 Tf-'^/^r veritical scale height of the gas disk 

V — acgh disk viscosity by a model 

/3 = — In Eg/ In a slope of the gas disk surface density 



Table 2: Statistics of the survival planets in 264 runs of simulation. 



Planet Mass 


No. 


Passing PAs 


Passing MMRs 


Final PAs 


Final MMRs 


Mp > 30Me 


473 


193 


2 


78 





lOMe < Mp < 30Me 


41 


18 


5 


12 





IMe < Mp < lOMe 


344 


297 


2 


148 





Mp < IMe 


579 


306 





82 





Total 


1437 


814 


9 


320 






Notes: PA means periastron alignment, MMR means mean motion resonance. Passing means 
at some stage of evolution (including the final time), planets are trapped either in PAs or 
MMRs. 
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Fig. 1. — Some tests are run to determine the appropriate values of type I migration reduction 
factor Ci in Eq.(l9]) and the eccentricity-damping rate K in Eq. (|T71) . (a) Migrations of a planet 
under different Ci . (b) Eccentricity evolutions of a planet under different K . 



-42 - 




Fig. 2. — Evolution of orbits in a typical run that two of the three survival planets passed 
through a periastron alignment. Initially 44 embryos are put with a disk depletion timescale 
Tdisk = 0.5 Myr. Three planets are survival with masses S.IM^, 238. 2M0, 137. GM^, and 
semimajor axes 1.07AU, 3.38AU, 6.98AU, respectively. (a)Planet growth and evolution in 
mass-semimajor axis plane. The red dots are the initial locations of the embryos. The green 
dots are those with either a > 50AU or e > 1 (being scattered out) during the evolution, 
(b) Evolution of semimajor axes of all embryos, (c) Evolution of periastron alignment angle 
{ix>2 — tJ^s). (d) Evolution of eccentricities for the three survival planets. 
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Fig. 3. — The influences of survival system properties due to different gas disks. The left 
vertical axis presents mean semimajor axis while the right one presents mean eccentricity, 
(a) influence of effective viscosity parameter a, (b) influence of the disk mass, the error bars 
are the standard deviations. 
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Fig. 4. — Distributions of semimajor axes for the survival planets in the 220 runs of sim- 
ulations at different epoches: (a) 0.01 Myr, (b) IMyr, (c) 3 Myr, (d) lOMyr during the 
evolution. 
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Fig. 5. — Correlations between the parameters of survival systems and the disk depletion 
timescale (xdisk) in 220 runs of simulations. Each point is average over 20 runs of systems with 
same r^isk, with error bars indicating the standard deviation. Panels from (a) to (d) shows 
the average eccentricities, numbers, masses and semimajor axes of the surviving planets 
versus Tdisk- 
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Fig. 6. — The correlations between the survival planet number A^ieft and the averaged ec- 
centricity Cmean (a) as Well as the average planet mass Mmean (b) in a planetary system. 
Both Cmean ciud Mmean are inversely correlated with A^ieft- The error bars are the standard 
deviations. 
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Fig. 7. — Distributions of semimajor axes, eccentricities and planet masses (a, e, Mp, resp.) 
from observations (left column, data: [http: / / exoplanet.eu ) and our simulations (right 
column). The shaded bars in (b),(d),(f) show distributions of the undetectable planets 
(with radial velocity Vr < 3m/s). The solid curves in (e) and (f) are fitting curves by 
A^(e)/iVtotai = 0.1v4exp(-y4e)/[l -exp(-y4)], with A = 4.2 for observations (e) and A = 7.8 
for simulations of planets with Vr > 3m/s (f). 
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Fig. 8. — Correlation graphes from observations(left column, data: http: / / exoplanet.euj) and 
our simulations (right column). Planets with induced stellar radial velocities Vr < 3m/s are 
shown in red triangles. Solid lines in (e) and (f) show the induced stellar radial velocities. 
The boxes in (c) & (d) show the planet desert in the Mp — e diagram. 
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Fig. 9. — Variations of mean eccentricities with semimajor axes (a) and distributions of 
eccentricities (b) for giant planets (Mp > lOM®) and terrestrial planets (Mp < lOM®). 
Histograms in (b) are fitted by A^(e)/A/'totai = 0.lAexp(-Ae)/[l - exp(-A)], with A = 7.78 
for giant planets and A = 21.0 for terrestrial planets. 
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Fig. 10. — Variations of eccentricity damping timescale [redap) for different planet masses 
(Mp) at semimajor axes 0.3AU, lAU, 3AU and 10 AU. For those planets Mp < Mj jj, T^dap 
is calculated by Eq.ffTOl) using a mean eccentricity e = 0.05. Otherwise r^dap is calculated 
by Eq. (fT7|) . The jump of T^dap occurs at Mp = Mjjj. Tedap keep horizontal until the mass 
of the planet becomes comparable to the disk mass. The Mp — e desert in Fig. 8d roughly 
corresponds to the mass regime with Tedap less than ~ IMyr. Out of this region eccentricities 
of planets can hardly be damped. 
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Fig. 11. — (a) The cumulative distribution function(CDF) of collision time, (b) The masses of 
embryos before collisions Mq and the masses after collisions Mi.(c) & (d): Mean semimajor 
axis and mass for TPs and GPs respectively. In panel (c): Mmean = 2.72M0, Omean = 
0.31AU (2D) and M^ean = 2.55Me,amean = 1.42AU (3D) for TPs. In panel (d): M_ = 
0.94Mj, a^nean = 3.41AU (2D) and M„ = 1.13Mj,a„ = 4.33AU (3D) for GPs. 
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Fig. 12. — Distributions of (a) semimajor axis, (b) planetary masses, (c) eccentricities 
and (d) inclinations of the survival planets in 20 runs of 2D simulations and 20 runs of 3D 
simulations, respectively. 
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Fig. 13. — Correlation graphes between inclinations with (a) semimajor axes, (b) masses in 
the 20 runs of 3D simulations, respectively. 



